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ABSTRACT 

Astronomical instruments generally possess spatially variant point-spread functions, which deter- 
mine the amount by which an image pixel is blurred as a function of position. Several techniques have 
been devised to handle this variability in the context of the standard image deconvolution problem. 
We have developed an iterative gravitational lens modeling code called Mirage that determines the 
parameters of pixelated source intensity distributions for a given lens model. We are able to include 
the effects of spatially variant point-spread functions using the iterative procedures in this lensing 
code. In this paper, we discuss the methods to include spatially variant blurring effects and test the 
results of the algorithm in the context of gravitational lens modeling problems. 
Subject headings: gravitational lensing: strong — methods: numerical 
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1. INTRODUCTION 

Modeling gravitational lens systems is fundamentally 
a two step process. A successful model requires an outer 
optimization procedure to discover the parameters of the 
lens model, and a nested optimization step to deter- 
mine the structure of the lensed source object. This in- 
ner source optimization step must incorporate the effect 
of the point spread function (PSF). Observed data are 
blurred by the PSF, and the presence of noise compli- 
cates the deconvolution process in general (|Hansen et al.l 
2006). This problem has been well studi ed by a variety 
of aut h ors using sp a tially i nvariant PSFs dWarren fc Dvel 
(p00l : IKoopmansI ISuvu et all (|2006D L 

To include blurring and lensing effects, we describe op- 
erations on images by the application of linear operators. 
We use "flattened" images to facilitate this notation, in 
which each column of the image is stacked upon the 
next, forming a vector. Our preferred scheme for solv- 
ing the lens modelin g problem is a modifi cation of the 
se milinear method of Warren fc Dvel (|2003D , as outlined 
in lRogers fc Fiegel (|2011[ ). To begin, we define image and 
source coordinate systems that are connected by the thin 
lens equation: 

P = 0-a(0), (1) 

where 9 and (3 are the image and source coordinates re- 
spectively, and a(9) is the deflection angle determined 
by the gravitational potential of the lens density dis- 
tribution. The lens equation is a nonlinear equation 
that maps pixel s from the image to the s ource plane 
([Schneider et al.l (fl992l ): iPetters et al.l (pOOl ). 

The source plane pixels are represented as a vector s, 
with elements that can be varied independently to ac- 
count for the details of the unknown source intensity dis- 
tribution. The details of the PSF are encoded in the blur- 
ring matrix B and gravitational lens effects, described 
by Equation [TJ are included in the lensing matrix L. We 
then define the total lens matrix / = BL and data vec- 
tor d. Denoting the standard deviation of the noise as 
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a, and given a source intensity distribution, the resulting 

x 2 



- 2 statistic between the data and model is 
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(2) 



By requiring derivatives of this equation with respect to 
the elements of the source vector s vanish, we find that 
the optimal source pixel intensities satisfy a least-squares 
equation: 

(3) 



F T Fs = F T d 



where we have absorbed factors of <Ti into the matrix 
Fij = fij/cri and data vector dj = dj I Oj . Th e deta ils of 
this derivat i on can be found in I Warren fc Dvel (|2003l ) and 
IKoopmansI (|2005f) . where the system is solved by direct 
matrix inversion with regularization. This least squares 
form is commonly found in the context of large scale 
image deconvolution pr oblems, which are t ypica ll y solved 
by ite r ative methods (jGolub fc Reinschl (|1970fl : [Biorck 
(fl996h :lH ansenl (|2010f) ) . Note that the semilinear method 
provides a technique for solving the linear parameters of 
the lensed system only (the source intensity distribution) 
in an "inner loop" , while the nonlinear parameters of the 
lens mass distribution m ust be solved in a separ ate "outer 
loop" optimization step (jRogers fc F iege 201l|). 

Spatial dependence of the PSF is not considered in 
most conventional deconvolution problems. This simpli- 
fies the construction of the blurring matrix B, since only 
one PSF is taken into account. However, it is well known 
that the PSF cannot always be treated as constant over 
an image in cases of astronomical interest. For example, 
spatially variant PSFs have be e n studied in the co ntext 
of adaptive optics (jLauerl (200i;[GnieseEaD(l2Q02)) and 
the PSF of astronomical instruments, such as the Hubble 
Advanced Camera for Surveys, ca n be extremely position 
dependent ()Bandara et al.l 12009). Several schemes have 
been designed to deal with this var i ability (IBoden et al.l 
([T995I) : iBirettal (|l99l : I Adorfl (|l99l : ILauerl (|2002fl ) . De- 
scribing a spatially variant PSF is much more compli- 
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cated than for the invariant case, since each row of the 
blurring matrix B will be derived from a unique PSF in 
general. The position of a pixel in the image determines 
the amount by which it is blurred. 

We illustrate the effect of spatially dependent PSFs on 
gravitationally lensed images in Figure [T] Consider the 
lensing effect produced by a singular isothermal sphere 
(SIS), which has three parameters: velocity dispersion 
a v and lens center (x,y). The deflection angle due to a 
SIS lens has a simple analytical form most conveniently 
described in standard polar coordinates (r,uj): 



a(r) =47r( — ) — - -, 
V c / D os 



(4) 



where D\ s and D os are the angular distances between lens 
and source and observer and source respectively, and c is 
the speed of light. We model the blurring in Figure[T]with 
a v = 265 km s^ 1 and use source redshift z s = 1.5 and 
lens plane redshift z\ = 0.12. This model was calculated 
using cosmological parameters H = 70 km s _1 Mpc -1 , 
Oo = 0.3 and Ao = 0.7 which we adopt for the remainder 
of this study. The source is comprised of a set of circular 
disks in the source plane as shown in the left-hand panel 
of Figure [1] and the gravitationally lensed image is shown 
in the center panel. The lensed image is then blurred by 
a spatially variant PSF and is shown in the right panel of 
Figure [TJ The distortion used to create this image varies 
from a delta function in the lower-left corner (negligible 
blur) to a Gaussian with standard deviation a g = 6.0 
pixels in the upper right corner. Each PSF is defined on 
an arbitrary 33 x 33 grid and is normalized to unity sum. 
The source and image plane size are 240 x 240 pixels. 

Unlike constant PSFs, spatially variant PSFs cannot 
be described by a simple convolution operation. For- 
tunately, numerical methods have been devised to han- 
dle th em, including sectioning methods (|Trussel fc Fogell 
1992), which deconvolve each PSF independently and 
forms the source from the sum of the reconstructions. 
iNaev fc O'Learvl (|1998f ) devised a clever method to 
model the effects of spatially variant PSFs within the 
framework of the standard image deconvolution problem. 
This approach differs from sectioning methods in that the 
separate PSFs are used to build an approximation to the 
blurred image of a given source, and a single iterative 
deconvolution operation is needed to solve for the source 
intensity di s tribut ion. The method was implemented in 
INaev et ail (|2002l ) and represents the spatial dependence 
of the PSF as a summation of piecewise blurring ma- 
trices, each of which applies over a limited area of the 
image . In this study, we use the method of INaev et al.l 
(2002) to incorporate spatially variant blurring into our 
gravitational lens modeling code. We briefly review the 
method here and discuss the procedure in detail in the 
Appendix. 

To include the effects of spatially variant blurs, the im- 
age of the unblurred l ensed source is pad ded to enforce 
a boundary condition (|Hansen et al.| [2006) . We focus on 
the use of reflexive boundary conditions, in which the im- 
age is padded by symmetric reflections of itself. Reflexive 
boundary conditions tend to reduce ringing artifacts if a 
significant amount of structure is located near the edges 



of the image. The image is then divided into a square 
grid, where the PSF is assumed constant in each region. 
These image regions and the PSFs are then padded to 
match in size. The two-dimensional fast Fourier trans- 
form (FFT) is used to calculate the resultant blurred im- 
age regions independently, resulting in an effective piece- 
wise convolution. By substituting efficient algorithms for 
the explicit matrix and matrix-transpose multiplications 
in Equation [3l the least squares form of the problem is 
preserved and the system can be solved efficiently. 

In principle spatially variant blurring can be described 
by a blurring matrix compatible with the semilinear 
method. However, in practice there are several prob- 
lems with the matrix approach. First, the size of the 
blurring matrix is N p i X x N p i X , so the matrix quickly 
becomes large as the image resolution is increased. Sec- 
ond, since the PSFs vary over regions of the image, it is 
possible that B may contain a large number of small but 
non-zero entries, particularly for large, complicated PSFs 
that are not well approximated by Gaussians or other 
simple analytical functions. This complicates the opti- 
mization because M—F T F must be inverted in the semi- 
linear scheme. It is generally required that M is sparse 
in order to store and invert this large matrix. The spar- 
sity requirement helps to reduce computation time and 
reduces the amplification of noise in the reconstructed 
source. In practice the semilinear method requires regu- 
larization to control the amount of noise present in the 
solution of Equation [3l The details and effects of several 
distinct regularization methods u sed with the semil inear 
method were studied i n detail bylSuvu et al l (|2006f ). 

Our previous work ((Rogers fc Fiegell201lD . compared 
the semilinear method with several iterative methods to 
solve the least-squares problem (Equation [3]). Iterative 
schemes have the advantage that time is saved by avoid- 
ing the explicit construction of the lens and blurring ma- 
trices. This is done using direct interpolation on the 
source plane under the effect of the lens equation (Equa- 
ti on HI). 



Roge rs fc Fiegd ( 20111 ) used the Qubist Optimization 
Toolbox (|Flegil20ll to map the x 2 surface over the 
space of the nonlinear lens parameters using the Ferret 
Genetic Algorithm (GA) and Locust Particle Swarm Op- 
timizer (PSO). Since this mapping requires a large num- 
ber of function evaluations (« 10 5 ) over the course of a 
run, speed is of the essence when choosing an inner loop 
optimization to determine the source plane parameters. 

Using the techniques introduced bv lNagv et al.1 (|2002l ) 
as a foundation, we have added the capability to include 
spatially variant PSFs to our gravitational lens modeling 
code using piecewise constant PSFs. This new capability 
is the subject of the current exploration. 

2. A SMALL-SCALE TEST 

In this section, we provide an example of modeling an 
extended source under the effects of a spatially variant 
PSF. We generate the lensed image of an analytical func- 
tion that describes a spiral source intensity distribution, 
according to the equation: 



5(r,w) = 



exp 



-2sin^ 



w 



rr 2 )] , (5) 
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where So is the maximum brightness in arbitrary units 
and core radius r c . The tightness of the arms about the 
central bulge is controlled by r, and ujq controls the orien- 
tation of the spiral, in standard polar coordinate s (r,ui). 
This a rtificial "galaxy" , originally described by IBonnet] 
(1995), serves as a convenient test pattern. To draw 
compa risons between our re sults for spatially invariant 
PSFs (|Rogers fe Fiegdl20ll . we w ill again make use of 
a Sing ular Isothermal Ellipse (SIE; iKeeton fc Kochanekl 
(1998)) with deflection angle components: 



bq 



: tan 



bq 



: tanh 



1>- 



q 2 s 



(6) 



(7) 



with V 2 = q 2 (s 2 + x 2 ) + y 2 and q = - e)/(l + e), 
and b is the equivalent Einstein radius when q = 1. The 
parameter b is related to the velocity dispersion a v by 
Equation HI 

The parameters used in this test are velocity dispersion 
a v = 265 km s~ 1 with Zd = 0.3 and z s = 1.05 giving an 
equivalent Einstein ring of b — 1.32 arcsec, ellipticity 
e = 0.35, lens center (x,y) = (0.11,0), core size s, and 
orientation angle 6l — 7r/4 measured counterclockwise 
from the right of the image. We set s = 0, resulting in a 
singular mass distribution. 

We used Equation [T] to form the lensed image of the 
source (Equation [5]) using the SIE deflection angle for- 
mulae. We generated a 20 x 20 grid of spatially variant 
Gaussian PSFs where each PSF is defined by a 33 x 33 
pixel mesh and has an FWHM ranging from 2.35 to 4.8 
pixels, shown in Figure [2] This grid of PSFs was used to 
blur the gravitationally lensed image and additive Gaus- 
sian white noise with standard deviation a g = 1.05 was 
added after the blurring operation, resulting in the arti- 
ficial data shown in Figure We define the peak signal- 
to-noise ratio (PSNR) as 



PSNR 



I„ 



(8) 



giving PSNR = 105.83. To illustrate the effect of vary- 
ing the number of PSFs, we model the data using smaller 
grids of 3 x 3, 5 x 5, 7 x 7, 10 x 10, and 20 x 20 PSFs. As 
shown in Figure |4j the best reconstruction with the low- 
est reduced \ 2 is obtained using a grid of 20 x 20 PSFs, 
which is the same number used to generate the data. 
This source and corresponding model image after 20 it- 
erations are also shown in Figure [3l The 3x3 and 5x5 
image residuals show significant structure, which is not 
present in the finer approximations. The residuals using 
a grid of 20 x 20 PSFs appear featureless. This demon- 
strates the improvement in image reconstruction as we 
include successively more information characterizing the 
blur. 

Figure [5] shows the relative error between the model 
source and the true solution as a function of iteration. 
For all PSF grid sizes, we find that the solutions display 
semi-convergence behavior such that the relative error 



between the model solution and the true solution im- 
proves until a minimum is reached and then begins to 
increase. This is due to the properties of the local opti- 
mizer used to determine the optimal source, and arises 
in the deconvolution step due to noise in the observed 
image. Regularization methods are generally used to 
control the increase of noise in t he reconstructed source 
found by the semilinear method (|Suvu et al.ll2006f ) . Sev- 
eral optimization methods have been applied to prob- 
lems with spatially variant blur i nc luding Landwebe r 
iteration (INocedal fc Wright! (fl999h ; iFish et all (fl99l : 
iTrussel fc Huntl (I1978D ). Richardson-Lucy deconvolution 
dFaisal et al.lll995|), and Lanczos-Tikhonov hybrid meth- 
ods (jChung et all [200*: in the context of the standard 
image deconvolution problem. Following Rog ers fc Fiegei 
(2011), we focus on the conjugate gradient method for 
least-squares problems (CGLS) and the steepest descent 
method (SD). Figure |5] shows the convergence history of 
th e SD algorithm. As in th e invariant PSF case discussed 
in iRogers fc Fiegei (|201 ID . the SD solution converges 
more slowly than CGLS and therefore it is less sensitive 
to the stopping criteria. When using an iterative method 
for local optimization, the number of iterations itself acts 
as a regularization parameter. The optimal stopping it- 
eration of these local optimizers is at the minimum of the 
relative error curve for a given set of lens parameters and 
PSF tiling. This critical iteration represents a balance 
between the re duced image y 2 an d the amount of regu- 
larization used ([Press et al.ll2007l ). Established methods 
exist to determine this critical iteration, i ncluding the 
L-Curve criterion ([Hansen fc Q'Learylll993f) and Gener- 
alized Cross Validation (IGolub et al.|[l979[) . In previous 
work (jRogers fc Fie gc 2011) we made use of the L-Curve 
criterion but Generalized Cross Validation is also imple- 
mented in our software. 

We find that the execution time of the problem in- 
cluding a spatially variant PSF increases approximately 
linearly with the number of separate PSFs used in the 
inversion as shown in Figure [71 This suggests that signif- 
icant gains could be made in the efficiency of the routine 
by parallelizing the implementation, since each image re- 
gion is independent. By splitting up the problem over 
several processors, the runtime for very large PSF grids 
can become feasible. 

3. A LARGE-SCALE TEST 

To demonstrate the code in operation on a large scale 
problem, we simulate the lensing effect of the mass dis- 
tribution of a galaxy cluster on a portion of the Hub- 
ble deep field using an elliptical potential. This test is 
intended as a demonstration of the feasibility and effi- 
ciency of our method on a problem that would be diffi- 
cult using the semilinear method while including a spa- 
tially dependent PSF. Problems of this size are realis- 
tic for a numbe r of pr actical modeling situations. For 
example, lAlardl ()2009[ ) has modeled the lensed system 
SL2SJ021408-053532, which produces a set of large arcs. 
This system has a lens that is comprised of a small group 
of six galaxies. Due to the large size of the lensed arcs, 
the scope of the source modeling prohibited the direct 
application of the semilinear method. 
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We form the lensed image of a portion of the Hubble 
deep field flWilliams et al.lll996l) by applying the ellipti- 
cal pote ntial oflBlandford &: Kochane k (1987) which was 
used bv lLink fc Piercd (|1998| ) to simulate the lens effect 
of the dark matter distribution of galaxy clusters. This 
potential function is given by 

52(1-9) 

y) = — ^— [s 2 + (1 + e c )x 2 + 2e s xy + (1 - e c )y 2 ] 9 , 

(9) 

which results in deflection angle at(0) — Vip(6). The 
elliptical potential depends on seven parameters: b is 
the equivalent Einstein radius in the limit of vanishing 
core radius s, ellipticity e, and power law index q, where 
< q < 0.5. The position angle of the lens <f> deter- 
mines the functions e c = ecos0 and e s = esin^>. We 
use the Einstein radius 6 = 9, power law index q = 0.25, 
<j> = it/A, position (x,y) — (0,0) and s — 0.5. The lens 
and source redshifts are Zd — 0.12 and z s = 1.5, re- 
spectively. We used an array of 25 PSFs arranged on 
a 5 x 5 grid to blur the image. This set of PSFs has 
been used to test image restoration schemes for Hub- 
ble Space Telescope (HST) images and represents the 
spatially variant nature of th e aberrations affect i ng th e 
HST before it was rep aired ([Katsaggcl os et al.l (|1994f ); 
iNagv k, O 'Learvl (|1998h ). The size of each PSF is 60 x 60 
pixels, and the source and image plane used to generate 
our lensed image are 800 x 800 pixel 2 . Gaussian white 
noise was added with standard deviation <j g = 1.37, giv- 
ing the image PSNR = 138.4. 

The image after 100 iterations is shown in Figure IH and 
a reduced x 2 = 0.995 was found. The system was solved 
using the CGLS algorithm with all 25 PSFs using the 
lens parameters defined above. The model took approx- 
imately 7 minutes to solve using a single 2.4 GHz CPU 
core. An approximation to the nonlinear lens parame- 



ters could be found using global optimization methods if 
one of the following strategies were employed: (1) a low- 
resolution approximation to the data could be used early 
during the lens parameter optimization, with successive 
refinement occurring later during the run; (2) a global 
optimizer could be used to roughly approximate the lens 
parameters, shifting to a faster local optimization scheme 
once solutions are localized to a small region of parame- 
ter space; or (3) global optimization could be used for the 
entire problem making use of large-scale parallelization. 

4. CONCLUSION 

We have developed a method to include the effects of a 
spatially variant PSF in gravitational lens modeling. In- 
cluding these effects in the standard semilinear method 
would be difficult due to the complicated blurring matrix 
required. These complicatio ns can be ov e rcome easily by 
incorporating the method of INagv et al.l (|2002t ) . Our ap- 
proach can acc ommodate la rge lensing problems like the 
case studied bv lAlardl (|2009f ). which limits the applicabil- 
ity of the direct semilinear approach. Techniques to in- 
clude the effects of spatially variant PSFs are important, 
as the response varies over the detector area for many 
astronomical instruments. Our algorithm allows this ef- 
fect to be included in lensing problems, thus improving 
the quality of reconstructions when the variability of the 
PSF is significant. The CGLS and SD algorithms allow 
a regularized inversion to be found quickly by truncated 
iteration. 
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APPENDIX 

SPATIALLY VARIANT BLURRING EFFECTS 

To describe blurring by a spatially variant PSF we first present an efficient method using two-dime nsional FFTs. We 
then s how how to treat the problem in term s of blurring matr ices and flattened image vectors. See INagv fc O 'Learvl 
(fl998h for more details on the approach and INagv et al.l (|2002h for a MATL AB implementation. 

Consider an N x N grid of independent PSFs and split the unknown blurred image Y into regions Yij , each of 
size k x k: 



Y u 


Y 12 




Yuf 


Y 21 


Y22 




Y~ 2N 










Y Nl 


Y N2 




Ynn 



Each of these blocks will be affected by an independent PSF. Suppose that the size of each PSF is (r + 1) x (r + 1) 
with r even, and let the unblurred N x N image be represented by X . 

Let us define a set of "mask" matrices Wij. In the case of piecewise constant PSFs, these masks are the same size 
as the unblurred image and are comprised of entries everywhere except for the k x k block at position where 
the entries of Wij are set to 1. 

To find the components of a given region we convolve X with the corresponding PSF , followed by an element- wise 
multiplication by the mask . The non-zero elements of this product give Y^ . Proceeding in this way we build up 
the blurred image block by block: 

N N 
»=1 3 = 1 
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where the symbol "o" represents element-wise multiplication and symbol "*" is the convolution operation. Note that 
each term in the sum is determined by the convolution of the entire image X with the appropriate PSF before the 
mask is applied. This is crucial to ensure that "seams" will not be visible between regions in the blurred image Y , 

In general, it is possible to speed up this routine by calculating Y^j directly. Consider splitting the unblurred image 
into regions X ^ where the superscript denotes the size of the block, in this case k x k. In order to avoid artifacts and 
keep the correct intensity near the edges of this block after convolution, we include a number of neighboring rows and 
columns on each side of X^. The width of this border is set by the size of the PSF, r/2, with regions on the image 
boundary padded to enforce the boundary conditions discussed in Section [T] These extended regions are then denoted 
X^ +k) . The PSFs are padded to match the extended regions in size, resulting in iy ■ The blurred extended region 
is found by the convolution 

rg +fc) = (pg +fc) *xg +fc) ). (A3) 

The central k x k block of this product is clipped out and placed in the position of Y. The process is repeated 
until the entire blurred image is filled in. Time is saved working with extended regions and padded PSFs since we 
only need to cal culate the convolution over the (r + k) x (r + k) block for each PSF rather than the entire image as 
in Equation IA21 and the construction of masks is not needed. The convolutions can be carried out efficiently with 
two-dimensional FFTs. 

The basic procedure can also be described by an analogous matrix-vector operation. To express the sum in Equation 

IA2l in terms of matrix multiplication, we define the unblurred flattened image as a vector x, and the flattened blurred 
image as y. We build a set of N 2 blurring matrices to describe the effect of each PSF on x, which we denote as By. 
The mask matrices Wij are used to construct analogous weighting matrices D ij . These matrices are of size N P i X x N P i X , 
where N p i X is the number of pixels in the image, identical to the size of the blurring matrices By. The total blurring 
matrix B is then written as a weighted sum of blurring matrices Bn, Bi2,...,Bjvi\r- 

N N 

B = EE D 'A- (A4) 
i=i j=i 

The blurred image is then found by a matrix multiplication y=Bx. The weighting matrices D, 3 have the mth diagonal 
entry equal to 1 provided that image pixel m is in region (i, j), and all other elements 0. The weighting matrices satisfy 

SjLi D%j = I where I is the N p i X x N p i x identity. We adopt the use of piecewise constant PSFs but in general 
it is possible to include higher order interpolation schemes between PSFs using the weighting matrices. Th e case of 
linear interpolation in solving systems with spatially variant blur has been studied by iNagy fc O'Learvl (|1998l) , but its 
inclusion complicates the pro cedure and did no t provide a significant improvement to the quality of the solution and 
increased computation times (jNagv et al.ll2002T ). 
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Source Plane 
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Blurred Image 




Fig. 1. — Example of spatially variant blurring. Left: a set of regular disks with radius 0.268 tile the source plane. Center: the circular 
disks are seen under the lensing effect of a Singular Isothermal Sphere (SIS) lens model. The SIS distorts the background circles into arcs, 
and the disk at the center of the SIS becomes a complete ring. Right: the same disk pattern under the effect of the SIS lens, with a spatially 
variant PSF blurring the observation. The blur is described by a delta function in the lower left hand corner to a Gaussian with standard 
deviation cr g = 6.0 pixels in the upper right corner, introducing a significant blur. 
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Fig. 2. — Grid of PSFs used in FigureQj] The PSFs vary from a Gaussian of FWHM of 2.35 pixels in the lower-left corner producing a 
modest blur to a Gaussian with FWHM 4.75 pixels in the upper right corner. 
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Fig. 3. — Top left: artificial data on a 120 X 120 grid. Bottom left: artificial source on a 50 X 50 grid. Top right: model observation. 
Bottom right: model source. The results after 19 iterations are shown. Note the presence of reconstructed noise in the source. The model 
has a reduced \ 2 = 0.998. 




Fig. 4. — Image residuals for a 3 X 3, 5 X 5, 7 X 7, 10 X 10, and 20 X 20 PSF grids after 19 CGLS iterations. For a small number of PSFs 
there is a significant amount of residual structure, but these artifacts are reduced as the grid of PSFs is enlarged. The reduced \ 2 is shown 
as a function of the number of PSFs (N pg f) used in the inversion. 
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Fig. 5. — Left: source convergence history using the CGLS algorithm. Right: corresponding Image convergence history. Note that the 
source displays semi-convergent behavior. The disagreement between model and actual source reaches a minimum before increasing. The 
critical iteration changes as the PSF grid is enlarged. 
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Fig. 6. — Left: source convergence history using the SD algorithm. Right: corresponding Image convergence history. The semi-convcrgent 
behavior of the source is less extreme than for the CGLS algorithm. 
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Fig. 7. — Timing results for the CGLS algorithm using N pa f as the number of PSFs to approximate the blurring effect. The plot illustrates 
the runtime for 4 X 4 to 20 X 20 square PSF grids in seconds. Each CGLS run was terminated at 20 iterations. 
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Fig. 8. — Top row: observation and model image. Middle row: actual and model source. The image and source plane are both 800 X 800 
pixels. These results are shown for 100 iterations. Bottom row: image residuals and an example of one of the 25 large PSFs used to generate 
the observations. Both of these images are plotted in logarithmic intensity to emphasize low level structure. Approximate runtime for this 
large-scale test is approximately 7 minutes. 



